#!/usr/bin/python
#Author: Jimmy Saw
#Date: 10-22-11

import sys
import os
import re
from Bio import SeqIO
from subprocess import call

def fmt(f):
    st = '{0:.4}'.format(f)
    return st

sctgfile = open(sys.argv[1], "rU")
sctgs = sctgfile.readlines()

mummerlist = []

for index, i in enumerate(sctgs):
    scaflist = []
    celera = open(i.strip()+".celerashreds.coords", "rU") #check Celera shreds
#    velvet = open(i.strip()+".velvetshreds.coords", "rU") #check Velvet shreds
    scaflist.append(i)
    cl = celera.readlines()
    tmp = cl[4]
    t = tmp.split('\t')
    ctgsize = int(t[7])
    ctgbegin = {}
    ctgend = {}
    for line in cl[4:]:
        l = line.split('\t')
        scafstart = int(l[0])
        scafstop = int(l[1])
        readname = l[-1].rstrip()
        scafname = l[-2]
        qcoverage = float(l[-3])
        if qcoverage < 100: #if alignment coverage is < 100, the rest is in another contig
            if scafstop < 1000: #work on beginning of contig
                ctgbegin[readname] = ((scafstart, scafstop, qcoverage)) #dictionary of tuples
            if scafstart > ctgsize - 1000: #work on end of contig
                ctgend[readname] = ((scafstart, scafstop, qcoverage))
    nctgbegin = {}
    nctgend = {}
    if index < len(sctgs)-1:
        nextscaf = sctgs[index+1]
        nc = open(nextscaf.strip()+".celerashreds.coords", "rU")
#        nv = open(nextscaf.strip()+".celerashreds.coords", "rU")
        ncl = nc.readlines()
        tmp2 = ncl[4]
        t2 = tmp2.split('\t')
        nctgsize = int(t2[7])
        for line in ncl[4:]:
            l = line.split('\t')
            nscafstart = int(l[0])
            nscafstop = int(l[1])
            nreadname = l[-1].rstrip()
            nscafname = l[-2]
            nqcoverage = float(l[-3])
            if nqcoverage < 100:
                if nscafstop < 1000:
                    nctgbegin[nreadname] = ((nscafstart, nscafstop, nqcoverage))
                if nscafstart > nctgsize - 1000:
                    nctgend[nreadname] = ((nscafstart, nscafstop, nqcoverage))

    for k, v in ctgend.iteritems():
        if k in nctgbegin:
            print "Between ", i.strip(), " and ", sctgs[index+1].strip(), k, fmt(v[2])            
#            tmpfasta = ""
#            call(["exfasta", k, "celera.shredded.ctgs.fasta"])

    celera.close()
    nc.close()
#generate list of gap-spanning shreds

#minimocmd = os.system("Minimo testgap_0069_0001.fasta -D QUAL_IN=testgap_0069_0001.qual\
#    -D MIN_LEN=30 -D FASTA_EXP=1 -D ACE_EXP=1 -D OUT_PREFIX=tmp_69_1")


